COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	\magnification\magstephalf
C00053 ENDMK
C⊗;
\magnification\magstephalf
\input macro.tex
\def\today{\ifcase\month\or
  January\or February\or March\or April\or May\or June\or
  July\or August\or September\or October\or November\or December\fi
  \space\number\day, \number\year}
\baselineskip 14pt
\rm
\line{\sevenrm a15.tex[162,phy] \today\hfill}

\parskip 5pt
\bigskip
\line{\bf Taylor\hfil}

\def\un{\underline}
\smallskip
Let $[q]$ be the characteristic function of a truth value~$q$;
$[F]=0$, $[T]=1$. An integral over a region, like 
$\int\int↓{(x,y)\in R}f(x,y)\,dy\,dx$, can be written as
$\int↓{x,y}[(x,y)\in R]f(x,y)$. An example application:
$$\displaylines{\int↓{xyz}[0<z<y<x<a]f(z)=\int↓z[0<z<a]\left(f(z)
\int↓{x,y}[z<y<x<a]\right)\cr
\noalign{\smallskip}
=\int↓z[0<z<a]f(z){(a-z)↑2\over 2}\,.\cr}$$

\proclaim Lemmas.
$$\displaylines{\int↓{x↓1,x↓2,\ldots,x↓n}[a<x↓1<x↓2<\cdots <x↓n<b]={(b-a)↑n\over
n!}\cr
\noalign{\smallskip}
\int↓{x↓0,x↓1,\ldots,x↓n}[a<x↓0<x↓1<\cdots <x↓n<b]f(x↓0)=\int↓{x↓0}
[a<x↓0<b]f(x↓0){(b-x↓0)↑n\over n!}\cr}$$
by the symmetry of the $n!$ possible orderings of $x↓1,x↓2,\ldots,x↓n$.

We next derive Taylor's theorem from the fundamental theorem of the calculus.
$$\eqalign{f(z)&=f(a)+\int↓y[a<y<z]f'(y)\cr
\noalign{\smallskip}
f'(y)&=f'(a)+\int↓x[a<x<y]f''(x)\cr
\noalign{\smallskip}
f''(x)&=f''(a)+\int↓w[a<w<x]f'''(w)\cr
\noalign{\smallskip}
f'''(w)&=f'''(a)+\int↓v[a<v<w]f''''(v)\,.\cr}$$
By making the obvious substitutions in the four instantiations of the
fundamental theorem, we get
$$\eqalign{f(z)&=f(a)+\int↓y[a<y<z]f'(a)+\int↓{xy}[a<y<z][a<x<y]f''(a)\cr
\noalign{\smallskip}
&\qquad\null+\int↓{wxy}[a<y<z][a<x<y][a<w<x]f'''(a)\cr
\noalign{\smallskip}
&\qquad\null+\int↓{vwxy}[a<y<z][a<x<y][a<w<x][a<v<w]f''''(v)\cr
\noalign{\smallskip}
&=f(a)+f'(a)\int↓y[a<y<z]+f''(a)\int↓{xy}[a<x<y<z]\cr
\noalign{\smallskip}
&\qquad\null+f'''(a)\int↓{xyz}[a<w<x<y<z]
+\int↓v[a<v<z]f''''(v)\int↓{wxy}[v<w<x<y<z]\cr
\noalign{\smallskip}
&=f(a)+(z-a)f'(a)+{(z-a)↑2\over 2!}f''(a)+{(z-a)↑3\over 3!}f'''(a)\cr
\noalign{\smallskip}
&\qquad\null+\int↓v[a<v<z]f''''(v){(z-v)↑3\over 3!}\,.\cr}$$
The last term is also equal to
$$\int↓{vwxy}[a<v<w<x<y<z]f''''(v)\,,$$
which for $z≥a$ is equal by the mean value theorem to
$$f''''(\zeta)\int↓{vwxy}[a<v<w<x<y<z]={(z-a)↑4\over 4!}f''''(\zeta)$$
for some $a<\zeta<z$.

\bigskip

\line{\bf Simpson\hfil}

\smallskip
The error $\epsilon$ in Simpson's rule is
$$\int↓z[-h<z<h]f(z)-{h\over 3}\bigl(f(-h)+4f(0)+f(h)\bigr)\,.$$
By the symmmetry around $z=0$, applying the rule to $g(z)={f(-z)+f(z)\over 2}$
gives the same error, so
$${\epsilon\over 2}=\int↓z[0<z<h]g(z)-{h\over 3}\bigl(2g(0)+g(h)\bigr)\,.$$
Plug in Taylor's theorem, noticing $g'(0)=g'''(0)=0$;
$$g(z)=g(0)+\int↓{xy}[0<x<y<z]g''(0)+\int↓{vwx}[0<v<w<x<y<z]g''''(v)$$
and cancel terms not involving $g''''$, because Simpson's rule is exact when
$g''''(v)\equiv 0$;
$$\eqalign{{\epsilon\over 2}&=\int↓{vwxyz}[0<v<w<x<y<z<h]g''''(v)
-{h\over 3}\int↓{vwxy}[0<v<w<x<y<h]g''''(v)\cr
\noalign{\smallskip}
&=\int↓v[0<v<h]g''''(v)\left({(h-v)↑4\over 4!}-{h\over 3}\,{(h-v)↑3\over 3!}\right)
\cr
\noalign{\smallskip}
&=-\int↓v[0<v<h]{(h-v)↑3(h+3v)\over 72}g''''(v)\cr
\noalign{\smallskip}
&=-\int↓v[0<v<h]w(v)g''''(v)\,,\quad w(v)≥0\,;\cr
\noalign{\smallskip}
{\epsilon\over 2}&=-g''''(\zeta)\int↓v[0<v<h]w(v)\cr
\noalign{\smallskip}
&=-g''''(\zeta)\left({h↑5\over 120}-{h\over 3}\left({h↑4\over 24}\right)\right)
={h↑5\over 180}g''''(\zeta)\,.\cr}$$
Subdivide the interval $(a,b)$ into $n$ parts and apply Simpson's rule to
each. Each part has width~$2h$, $h={b-a\over 2n}$; the total error is
$${h↑5\over 360}\sum g''''(\zeta↓i)={h↑4\over 720}\sum 2hg''''(\zeta↓i)\,,$$
where the sum approaches $\int↓n[a<u<b]g''''(u)=g'''(b)-g'''(a)$; the
error is therefore asymptotically
${(b-a)↑4\over 5760n↑3}\bigl(g'''(b)-g'''(a)\bigr)$.

Simpson's rule approximates an integral by a sum, and estimates the error.
One can equally approximate a sum by an integral and estimate the error.

\bigskip
\line{\bf Euler (I)\hfil}

\smallskip
The difference 
$$D=\int↓z[0<x<1]f(z)-{f(0)+f(1)\over 2}$$
is zero if $f$ is a polynomial of degree~1, i.e., if $f''≡0$. This and
the linearity of~$D$ suggest that $D$ can be expressed as a weighted
integral of~$f''$. To find the formula, replace~$f$ in the formula
for~$D$ by its Taylor expansion with exact remainder in~$f''$.
$$f(z)=f(0)+\int↓y[0<y<z]f'(0)+\int↓{x,y}[0<x<y<z]f''(x)\,.$$
Then
$$\eqalign{D&=\int↓z[0<z<1]f(0)+\int↓{yz}[0<y<z<1]f'(0)+
\int↓{xyz}[0<x<y<z<1]f''(x)\cr
\noalign{\smallskip}
&\qquad\null-{f(0)\over 2}\cr
\noalign{\smallskip}
&\qquad\null-{1\over 2}\,f(0)-{1\over 2}\int↓y[0<y<1]f'(0)-{1\over 2}
\int↓{xy}[0<x<y<1]f''(x)\,.\cr}$$
The terms not containing $f''$ vanish, whether by manipulation or by the
above argument. Then
$$\eqalignno{D&=\int↓{xyz}[0<x<y<z]f''(x)-{1\over 2}\int↓{xy}[0<x<y<1]f''(x)&*1\cr
\noalign{\smallskip}
&=\int↓x[0<x<1]
\left({(1-x)↑2\over 2}-{(1-x)\over 2}\right)\,f''(x)\cr
\noalign{\smallskip}
&=\int↓x[0<x<1]
{x(x-1)\over 2}\,f''(x)\,.&*2\cr}$$
Because the weight ${x(x-1)\over 2}$ does not change sign in the region of 
integration, the mean-value theorem gives
$$D=f''(\zeta)\int↓x[0<x<1]{x(x-1)\over 2}=-{f''(\zeta)\over 12}\qquad
\hbox{for some $0<\zeta <1$}.\eqno{*3}$$
If $f↑{(4)}≡0$, $D$ is $-{1\over 12}$ times the average of~$f''$, or
$-{1\over 12}\bigl(f'(1)-f'(0)\bigr)$. This suggests expressing~$D$
as $-{1\over 12}\bigl(f'(1)-f'(0)\bigr)$ plus a weighted integral of~$f↑{(4)}$.
Using *1, we substitute
$$f''(x)=f''(0)+\int↓w[0<w<x]f'''(0)+\int↓{vw}[0<v<w<x]f↑{(4)}(v)$$
and
$$f'(y)=f'(0)+\int↓x[0<x<y]f''(0)+\int↓{wx}[0<w<x<y]f'''(0)+\int↓{vwx}
[0<v<w<x<y]f↑{(4)}(v)$$
into
$$\eqalign{D↓2&=D+{1\over 12}\bigl(f'(1)-f'(0)\bigr)\cr
\noalign{\smallskip}
&=\int↓{xyz}[0<x<y<z<1]f''(x)-{1\over 2}\int↓{xy}[0<x<y<1]f''(x)+{1\over 12}
\bigl(f(1)-f(0)\bigr)\,,\cr}$$
yielding
$$\eqalign{D↓2&=\int↓{xyz}[0<x<y<z<1]f''(0)+\int↓{wxyz}[0<w<x<y<z<1]f'''(0)\cr
\noalign{\smallskip}
&\qquad\null+\int↓{vwxyz}[0<v<w<x<y<z<1]f↑{(4)}(v)-{1\over 2}\int↓{xy}
[0<x<y<1]f''(0)\cr
\noalign{\smallskip}
&\qquad\null-{1\over 2}\int↓{wxy}[0<w<x<y<1]f'''(0)-{1\over 2}\int↓{vwxy}
[0<v<w<x<y<1]f↑{(4)}(v)\cr
\noalign{\smallskip}
&\qquad\null+{1\over 12} f'(0)+{1\over 12}\int↓x[0<x<1]f''(0)+{1\over 12}
\int↓{wx}[0<w<x<1]f'''(0)\cr
\noalign{\smallskip}
&\qquad\null+{1\over 12}\int↓{vwx}[0<v<w<x<1]f↑{(4)}(v)-{1\over 12}f'(0)\cr}$$
where all terms not containing $f↑{(4)}$ cancel;
$$\eqalignno{D↓2&=\int↓{vwxyz}[0<v<w<x<y<z<1]f↑{(4)}(v)-{1\over 2}
\int↓{vwy}[0<v<w<x<y<1]f↑{(4)}(v)\cr
\noalign{\smallskip}
&\qquad\null+{1\over 12}\int↓{vwx}[0<v<w<x<1]f↑{(4)}\cr
\noalign{\smallskip}
&=\int↓v[0<v<1]\left({(1-v)↑4\over 24}-{1\over 2}\,{(1-v)↑3\over 6}+
{1\over 12}\,{(1-v)↑2\over 2}\right)f↑{(4)}(v)\cr
\noalign{\smallskip}
&=\int↓v[0<v<1]f↑{(4)}(v){v↑2(1-v)↑2\over 24}&*4\cr
\noalign{\smallskip}
&=f↑{(4)}(\zeta)/720\,,\qquad\hbox{for some $0<\zeta<1$}.\cr}$$
Putting results together,
$$\int↓z[0<z<1]
f(z)-{f(0)+f(1)\over 2}=-{1\over 12}\,\bigl(f'(1)-f'(0)\bigr)
+\int↓v[0<v<1]f↑{(4)}(v)\,{v↑2(1-v)↑2\over 24}\,.$$
Applying analogous formulas to the intervals $(1,2),(2,3),\ldots,(n-1,n)$,
and adding
$$\eqalign{%
\int↓z [0<z<n]
f(z)&-\left({f(0)\over 2}+f(1)+f(2)+\cdots+f(n-1)+{f(n)\over 2}
\right)\cr
\noalign{\smallskip}
&=-{1\over 12}\bigl(f'(n)-f'(0)\bigr)+\int↓v[0<v<1]f↑{(4)}(v)w(v)\,.\cr}$$
If either the integral or the sum on the left of the equation is known,
the other can be estimated. Use this to estimate 
$\ln(n!)=\ln(1)+\ln(2)+\cdots+\ln(n)$; 
$$\int↓z[0<z<n]\ln(z)-\ln(n!)+{\ln n\over 2}=-{1\over 12}\left({1\over n}
-1\right)+\int↓v[0<v<n]{w(v)\over v↑4}\,,$$
where the first integral is $n\ln n-n+1$ and the second is $O(1)$.
$$\eqalign{\ln(n!)&=n\ln n-n+{\ln n\over 2}+O(1)\cr
\noalign{\smallskip}
n!&=\left({n\over e}\right)↑n\sqrt{n}\,O(1)\,,\cr}$$
a simple form of Stirling's approximation. The above derivation is a special
case of Euler's summation formula.

\medskip
\line{\bf Uspenskii\hfil}

\smallskip
Next we estimate the central binomial coefficient ${2n\choose n}$ for large~$n$,
using only elementary means, following Dynkin and Uspenskii.

Let $x↓n$ be the central term in the binomial distribution for $p=1/2$:
$$\eqalign{x↓n&={2n\choose n}\,2↑{-2n}={1\cdot 3\cdot 5\cdot \;\cdots \;\cdot
(2n-1)\over 2\cdot 4\cdot 6\cdot\; \cdots\; \cdot (2n)}\cr
\noalign{\smallskip}
x↓0&=1\,;\quad x↓1=1/2\,;\quad x↓2=3/8\,;\quad x↓4=35/128\,.\cr}$$
Clearly $x↓{n+1}={2n+1\over 2n+2}\,x↓n$, $x↓n\searrow 0$ (decreases to~0 
as limit). We shall see that $x↓n\in\Theta(n↑{-1/2})$; in fact,
$\pi nx↓n↑2→1$. Define $y↓n=nx↓n↑2<v↓n=(n+1/4)x↓n↑2<u↓n=
{n↑2\over n-1/4}\,x↓n↑2<z↓n={n↑2\over n-1}\,x↓n↑2$. Because of the
recurrence relations
$$\eqalign{{y↓{n+1}\over y↓n}&={(2n+1)↑2\over 2n(2n+2)}=1+{1\over 4n↑2+4n}\cr
\noalign{\smallskip}
{v↓{n+1}\over v↓n}&={(4n+5)(4n+2)↑2\over (4n+1)(4n+4)↑2}=1+{1\over 16n↑3+36n↑2+24n+4}\cr
\noalign{\smallskip}
{u↓{n+1}\over u↓n}&={(4n-1)(4n+2)↑2\over (4n+3)(4n)↑2}=1-{1\over 16n↑3+12n↑2}\cr
\noalign{\smallskip}
{z↓{n+1}\over z↓n}&={(2n-2)(2n+1)↑2\over (2n)↑3}=1-{3n+1\over 4n↑3}\cr}$$
we can see that $y↓n$ and $v↓n$ increase, $u↓n$ and $z↓n$ decrease. Because
$z↓n/y↓n=n/(n-1)\searrow 1$, all four have a common limit~$L$. It is not hard
to see by the recurrence relations that $(3y↓n+z↓n)/4$ converges faster than
$y↓n$ or~$z↓n$, and that $(v↓n+u↓n)/2$ converges faster than $v↓n$ or~$u↓n$.

\bigskip
$$\vcenter{\halign{$\hfil#\hfil$
&$\hfil#\hfil$\qquad
&$\hfil#\hfil$
&$\hfil#\hfil$\quad
&$\hfil#\hfil$
&$\hfil#\hfil$\quad
&$\hfil#\hfil$
&$\hfil#\hfil$\qquad
&$\hfil#\hfil$
&$\hfil#\hfil$\cr
y↓n&&v↓n&&&&&u↓n&&z↓n\cr
\bullet&&\bullet&&\bullet&&&\bullet&&\bullet\cr
&&&&&L\cr
&y↓{n+1}&&v↓{n+1}&&&u↓{n-1}&&z↓{n-1}\cr
&\bullet&&\bullet&\bullet&&\bullet&&\bullet\cr
}}$$

\bigskip
For $a≤b$, the monotonicity of $y↓n$, $v↓n$, $u↓n$, and $z↓n$ implies:
$$\eqalign{\sqrt{ax↓a↑2}&≤\sqrt{b}\,x↓b≤\sqrt{L}\cr
\noalign{\smallskip}
\sqrt{(a+1/4)x↓a↑2}&≤\sqrt{b+1/4}\,x↓b≤\sqrt{L}\cr
\noalign{\smallskip}
\sqrt{L}&≤\sqrt{b↑2/(b-1/4)}\,x↓b≤\sqrt{(a↑2x↓a↑2)/(a-1/4)}\cr
\noalign{\smallskip}
\sqrt{L}&≤\sqrt{b↑2/(b-1)}\,x↓b≤\sqrt{(a↑2x↓a↑2)/(a-1)}\,.\cr}$$
Solving the second and third inequalities for $x↓b$ in terms of~$L$, 
the resulting bounds
$$\sqrt{L(b-1/4)/b↑2}≤x↓b≤\sqrt{L/(b+1/4)}$$ 
are in the ratio $\sqrt{(b↑2-1/16)/b↑2}$, which is less than~1 but greater
than $1-{1\over 32b↑2}$. A~geometric or arithmetic mean of the two bounds is
therefore a quite accurate estimate of~$x↓b$.

Using the data for $a=34$ in the table below shows that for $b≥34$,
$$\displaylines{\sqrt{0.3159780400/b}≤x↓b\cr
\noalign{\smallskip}
\sqrt{0.3183014080/(b+1/4)}≤x↓b\cr
\noalign{\smallskip}
x↓b≤\sqrt{0.3183186181(b-1/4)/b↑2}\cr
\noalign{\smallskip}
x↓b≤\sqrt{0.3255531321(b-1)/b↑2}\,.\cr}$$

{}
$$\vcenter{\offinterlineskip
\halign{\quad$\hfil#\hfil$\quad%
&\vrule#&\strut\quad$#\hfil$\quad%
&\vrule#&\quad$#\hfil$\quad%
&\vrule#&\quad$#\hfil$\quad\cr
a&&\hfil 10&&\hfil 25&&\hfil 34\cr
&height2pt&\omit&&\omit&&\omit\cr
\noalign{\hrule}
&height2pt&\omit&&\omit&&\omit\cr
x↓a&&0.1761970520&&0.1122751727&&0.09640265435\cr
&height2pt&\omit&&\omit&&\omit\cr
x↓a↑2&&0.03104540113&&0.0126057144&&0.009293471765\cr
&height2pt&\omit&&\omit&&\omit\cr
y↓a&&\un{0.31}04540113&&\un{0.31}51428600&&\un{0.31}59780400\cr
&height2pt&\omit&&\omit&&\omit\cr
v↓a&&\un{0.3182}153616&&\un{0.31829}42885&&\un{0.31830}14080\cr
&height2pt&\omit&&\omit&&\omit\cr
u↓a&&\un{0.3184}143706&&\un{0.3183}261212&&\un{0.31831}86181\cr
&height2pt&\omit&&\omit&&\omit\cr
z↓a&&\un{0.3}449489015&&\un{0.32}82738124&&\un{0.32}55531321\cr
&height2pt&\omit&&\omit&&\omit\cr
{v↓a+u↓a\over 2}&&\un{0.31831}48661&&\un{0.318310}2048&&\un{0.318310}0131\cr
&height2pt&\omit&&\omit&&\omit\cr
{3y↓a+z↓a\over 4}&&\un{0.319}0777339&&\un{0.3184}255981&&\un{0.3183}718130\cr
&height2pt&\omit&&\omit&&\omit\cr
}}$$

\noindent
Since $1/\pi=0.3183098862$, the value of $(v↓{34}+u↓{34})/2$ suggests
$L=1/\pi$. Underlining indicates number of digits of agreement with~$L$.

To find the limit $L$, and other useful stuff, expand
$$\eqalignno{L&=\lim {2n\over 2}\,x↓n↑2=\lim {1\over 2}\times{1\over 2}
\times {1\over 2}\times{3\over 4}
\times {3\over 4}\times{5\over 6}
\times {5\over 6}\times\cdots\times{2n-1\over 2n}\times (2n-1)\cr
\noalign{\smallskip}
&={1\over 2}\times{1\over 2}
\times {3\over 2}\times{3\over 4}
\times {5\over 4}\times{5\over 6}
\times {7\over 6}\times{7\over 8}\times
\cdots\times{2n-1\over 2n-2}\times {2n-1\over 2n}\times\cdots&(*)\cr
\noalign{\smallskip}
&=1/\pi\qquad\hbox{(Wallis, 1655)}\cr
\noalign{\smallskip}
&={1\over 4}\times{9\over 8}
\times {25\over 24}\times{49\over 48}\times\cdots\cr
\noalign{\smallskip}
&={1\over 2}\times{3\over 4}
\times {15\over 16}\times{35\over 36}\times{63\over 64}\times\cdots\;.\cr}$$
The starred line may be treated as the product analogue of an alternating
series. For example, if $S↓i$ is the $i$-th partial product,
$S'↓i\defeq\sqrt{S↓iS↓{i+1}}$ has the same limit, converges faster, and
can be treated again the same way. Using $S↓{10}$ through $S↓{16}$ and
transforming the series six times gives an estimate 0.3183098559; the
reciprocal is $3.141592953=\pi$ to~7S.

Example: Flip a fair coin two million times. The chance that it comes
up heads exactly half the time is bounded by
$$\sqrt{{10↑6-.25\over 10↑{12}\pi}}≤x↓{10↑6}≤\sqrt{{1\over\pi(10↑6+.25)}}\,.$$
Both bounds, to 10S, are $5.641895130\times 10↑{-4}$. 

Flip it a thousand times:
$$\displaylines{\sqrt{{500-.25\over 250000\pi}}≤x↓{500}≤(500.25\pi)↑{-1/2}\cr
\noalign{\smallskip}
0.02522501660≤x↓{500}≤0.02522501975\cr}$$
accurate to~7S. Typically with asymptotic expansions, the larger~$n$ is,
the easier it is to calculate~$x↓n$ accurately. In this problem,
cases around~$x↓{100}$ are the hardest to evaluate to high precision.

It is known (by Euler's summation formula, demonstrated below) that
$n!\,e↑n\,n↑{-(n+1/2)}$ has a limit,~$K$. Then
$$y↓n=n\,x↓n↑2={n(2n)!↑2\over n!↑4\,2↑{4n}}=
{\bigl((2n)!\,e↑{2n}\,(2n)↑{-(2n+1/2)}\bigr)↑2\over
(n!\,e↑n\,n↑{-(n+1/2)})↑4}\cdot 2→{2\over K↑2}$$
and $y↓n→L=1/\pi$, so $K=\sqrt{2\pi}$. Therefore
$$n!=\left({n\over e}\right)↑n\,\sqrt{{n\over 2\pi}}\,\bigl(1+o(1)\bigr)\,.$$

What has gone before is intended, by using improvised methods on special
cases, to suggest an approach to approximating large sums closely by
integrals, resulting in Euler's summation formula.

\vfill\eject

\line{\bf Euler $+$ Bernoulli $=$ Stirling $+$ Glaisher\hfil}

Suppose $p↓j(x)$ is an integrable function of period~1. In general,
$\int↓0↑xp↓j(y)\,dy$ is not a periodic function of~$x$, but if we subtract
from~$p↓j$ its average value $m↓j\defeq\int↓0↑1p↓j(y)\,dy=\int↓x↑{x+1}p↓j(y)\,dy$,
we get $q↓j(x)\defeq p↓j(x)-m↓j$, where $q↓j(x)$ is periodic and so is its
negative integral $p↓{j+1}(x)\defeq -\int↓0↑xq↓j(y)\,dy$. Starting with any
integrable function as~$p↓1$, induction defines the two infinite families
$p↓j(x)$ and $q↓j(x)$, $1≤j$. The integrations make them successively
smoother: $p↓j$~and $q↓j$ have at least $j-2$ continuous derivatives. They
satisfy the equations
$$\eqalign{p↓{j+1}(0)&=0\cr
\noalign{\smallskip}
p'↓{j+1}(x)&=q'↓{j+1}(x)=-q↓j(x)\cr
\noalign{\smallskip}
\int↓x↑{x+1}q↓j(y)\,dy&=0\,.\cr}$$

Specifying $p↓1$ in more detail determines more properties of families~$\{p↓j\}$
and~$\{q↓j\}$.

\display 20pt:$\bullet$:
If $p↓j$ is a polynomial of degree~$d$ in the interval $U=[0,1)$, so is~$q↓j$,
while $p↓{j+1}$ is a polynomial of degree $d+1$ in~$U$.

\display 20pt:$\bullet$:
If $p↓1$ is sinusoidal, say $p↓1=\exp(2\pi kxi)$, $i=\sqrt{-1}$, then all the
$p↓j$ and~$q↓j$ are sinusoidal with the same period $1/k$, with amplitude 
proportional to $(2\pi k)↑{-j}$. For any nonconstant~$p↓1$, for large~$j$, 
$q↓j$~is very nearly sinusoidal.

\display 20pt:$\bullet$:
If $p↓j$ is even, it is symmetric in $U$: $p↓j(x)=p↓j(-x)=p↓j(1-x)$. Similarly,
if $p↓j$ is odd, it is antisymmetric in~$U$:
$p↓j(x)=-p↓j(-x)=-p↓j(1-x)$. If $p↓j$ is odd, $m↓j=0$, so if $p↓j$ is
odd/even, $q↓j$~is odd/even and $p↓{j+1}$ is even/odd.

\display 20pt:$\bullet$:
For $j>1$, $p↓j(0)=p↓j(1)=0$. For convenience and uniformity, without
changing any of the other functions we may add a
constant to~$p↓1$ so that $p↓1(0)=p↓1(1)=0$ also.

\display 20pt:$\bullet$:
Suppose $p↓1$ is odd, with finitely many discontinuities in~$U$, and is strictly
positive in $L=(0,1/2)$, negative in $R=(1/2,0)$.

\vfill

\centerline{Typical $p↓1=q↓1$}

\eject

\display 20pt::
Then $m↓1=0$, $q↓1=p↓1$; being the integral of~$q↓1$, $p↓2$~is even and strictly
increasing in~$L$, strictly decreasing in~$R$, with minima at $(0,0)$ and
$(1,0)$, a~unique maximum in~$U$ at $x=1/2$, and strictly positive in $(0,1)$;
$m↓2>0$.

\figbox 2.5in:

\centerline{Typical $p↓2$\qquad%
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad And $q↓2$}

\display 20pt::
Then $q↓2$ is even, strictly increasing/decreasing in $L/R$, has a unique
and negative minimum in~$U$ at $x=0$, a~unique and positive maximum at
$x=1/2$. It has two roots, $r\in L$ and $1-r\in R$. The half-integrals
$\int↓0↑{1/2}q↓2(y)\,dy$ and $\int↓{1/2}↑1 q↓2(y)\,dy$ are both zero.

\figbox 2.5in:

\centerline{Typical $p↓3=q↓3$}

\medskip
\display 20pt::
Then $p↓3=q↓3$ is odd, strictly decreasing in $(0,r)$, increasing in $(r,1-r)$,
and decreasing in $(1-r,1)$; has a unique and negative minimum in~$U$ at~$r$,
a~unique and positive maximum at $1-r$, with roots in~$U$ only at~0 and $1/2$.
Because  $-p↓3$ meets the same conditions as~$p↓1$, an obvious induction gives
similar properties for all $\{p↓j\}$ and $\{q↓j\}$.

\vfill\eject

When we integrate a function $f$ between integer limits $(a,b)$ with
weight~$q↓j$, we find
$$\eqalign{\int↓a↑bf(x)q↓j(x)\,dx&=-\int↓a↑bf(x)p'↓{j+1}(x)\,dx\cr
\noalign{\smallskip}
&=-\bigl(f(x)p↓{j+1}(x)\bigr)↓a↑b+\int↓a↑bf'(x)p↓{j+1}(x)\,dx\cr
\noalign{\smallskip}
&=0+m↓{j+1}\int↓a↑bf'(x)\,dx+\int↓a↑bf'(x)q↓{j+1}(x)\,dx\cr
&=m↓{j+1}\bigl(f(x)\bigr)↑b↓a+\int↓a↑bf'(x)q↓{j+1}(x)\,dx\,.\cr}$$
Inductively,
$$\int↓a↑bf(x)q↓j(x)\,dx=\sum↓{0≤\lambda <k}m↓{j+\lambda +1}
\bigl(f↑{(\lambda)}(x)\bigr)↓a↑b+\int↓a↑b
f↑{(k)}(x)q↓{j+k}(x)\,dx\,.\eqno(*)$$

What follows examines the simplest nontrivial choice for~$p↓1$, where
$p↓1(x)=x-\lfloor x\rfloor$, the fractional part of~$x$; on~$U$,
$p↓1(x)=x$. Tabulating, (see graphs),
$$\vcenter{\halign{$\hfil#\hfil$\qquad
&$#\hfil$\qquad&$\hfil#\hfil$\qquad
&\hfil#\cr
&&&\hbox{lcm of denomi-}\cr
j&p↓j\hbox{ (in $U$)}&m↓j&\hbox{nators in $q↓j$}\hfil\cr
\noalign{\medskip}
1&x&1/2&2\cr
\noalign{\smallskip}
2&-{x↑2\over 2}+{x\over 2}&1/12&12\cr
\noalign{\smallskip}
3&{x↑3\over 6}-{x↑2\over 4}+{x\over 12}&0&12\cr
\noalign{\smallskip}
4&-{x↑4\over 24}+{x↑3\over 12}-{x↑2\over 24}&-1/720&24\cr
\noalign{\smallskip}
5&{x↑5\over 120}-{x↑4\over 48}+{x↑3\over 72}-{x\over 720}&0&720\cr
\noalign{\smallskip}
6&-{x↑6\over 720}+{x↑5\over 240}-{x↑4\over 288}+{x↑2\over 1440}&1/30240&1440\cr
\noalign{\smallskip}
7&{x↑7\over 5040}-{x↑6\over 1440}+{x↑5\over 1440}-{x↑3\over 4320}+{x\over 30240}%
&0&30240\cr
\noalign{\smallskip}
8&-{x↑8\over 40320}+{x↑7\over 10080}-{x↑6\over 8640}+{x↑4\over 17280}-{x↑2\over
60480}&-1/1209600&120960\cr
}}$$
These are the related to the Bernoulli numbers $B↓j$ and polynomials
$B↓j(x)$:
$$\eqalign{m↓j&=(-)↑jB↓j/j!\,,\qquad\hbox{where the sign affects only $m↓1$,}\cr
\noalign{\smallskip}
q↓j(x)&=(-)↑jB↓j(x)/j!\,.\cr}$$

\vfill\eject

figure goes here

\vfill\eject

To accompany this choice of $\{p↓j,q↓j\}$ for $j≥1$, what $p↓0$ and $q↓0$
would obey the same laws? We would have $q↓0(x)=-p'↓1(x)=-1$ for $0<x<1$, but
at $x=1$ the integral of~$q↓0$ must go discontinuously from $-1$ to~0.
This requires $q↓0(x)$ to have, in effect, a~sharp peak just below $x=1$,
with negligible width and unit area; $p↓0(x)$ is greater by~1, to have
$p↓0(0)=0$, so $m↓0=1$.

\bigskip
\bigskip
$$\vcenter{\halign{$\hfil#\hfil$\qquad%
&$\hfil#\hfil$\qquad%
&$\hfil#\hfil$\enspace%
&$\hfil#\hfil$\qquad%
&$\hfil#\hfil$\enspace%
&$\hfil#\hfil$\quad%
&$\hfil#\hfil$\quad%
&$\hfil#\hfil$\enspace%
&$\hfil#\hfil$\enspace%
&$\hfil#\hfil$\cr
&&&&&0&&&&1\cr
&&.&\bullet&0&.&&.&.\cr
0&&&1&-1&.&&.&.\cr
&p↓0&&&&&&q↓0\cr
}}$$
These fictional functions, approximated by the graphs, have the very useful
property of relating sums to integrals:
$$\eqalign{\int↓a↑bf(x)p↓0(x)\,dx&=\sum↓{a<i≤b}f(i)\cr
\noalign{\smallskip}
\int↓a↑bf(x)q↓0(x)\,dx&=\sum↓{a<i≤b}f(i)-\int↓a↑bf(x)\,dx\,.\cr}$$

Conveniently forgetting that there are no such functions, we ``deduce''
$$\eqalignno{\sum↓{a<i≤b}f(i)-\int↓a↑bf(x)\,dx&=\int↓a↑bf(x)q↓0(x)\,dx\cr
\noalign{\smallskip}
&=m↓1\bigl(f(x)\bigr)↓a↑b+\int↓a↑bf'(x)q↓1(x)\,dx\,.&(**)\cr}$$
The first formula is, in fact, equal to the last, by induction on~$b$. The
induction step shows that
$$f(b+1)-\int↓b↑{b+1}f(x)\,dx={f(b+1)-f(b)\over 2}+\int↓b↑{b+1}f'(x)(x-b-1/2)\,dx$$
by partial integration of the second integral.

Take some fixed function $f$, with $i$-th derivative $f↑{(i)}$ ($1≤i≤k$). Define
$f↑{(0)}(x)=f(x)$, $f↑{(-1)}(x)=$ any antiderivative (indefinite integral)
of $f(x)$. Let $F(x)$ be a summation of $f(x)$, i.e., $F(x)=F(x-1)+f(x)$. Define
$$\sigma↓{ab}=\bigl(F(x)\bigr)↓a↑b=\sum↓{a<i≤b}f(i)\,.$$
Then, combining $(*)$ and $(**)$ gives

\medskip\noindent
{\it Euler's summation formula:\/} $\sigma↓{ab}=\bigl(g↓k(x)\bigr)↓a↑b+E↓k(a,b)$

\smallskip\noindent
where the approximate indefinite summation is
$$\eqalign{g↓k(x)&=\sum↓{0≤j≤k}m↓jf↑{(j-1)}(x)=\sum↓{0≤j≤k}\,T↓j\cr
\noalign{\smallskip}
&=f↑{(-1)}(x)+{f(x)\over 2}+{f'(x)\over 12}-{f'''(x)\over 720}+{f↑{(5)}\over
30240}+\cdots+m↓kf↑{(k-1)}(x)\cr}$$
and the exact error of the approximation is
$$E↓k(a,b)=\int↓a↑bf↑{(k)}(x)q↓k(x)\,dx\,.$$

If $k=1$, Euler's formula gives
$$f(a+1)+f(a+2)+\cdots +f(b)=\int↓a↑bf(x)\,dx+{1\over 2}\bigl(f(b)-f(a)\bigr)
+\int↓a↑bf'(x)(\{x\}-1/2)\,dx$$
which is the trapezoidal rule with exact error term:
$$\int↓a↑bf(x)\,dx={f(a)\over 2}+f(a+1)+f(a+2)+\cdots+f(b-1)+{f(b)\over 2}
+\int↓a↑bf'(x)(1/2 \{x\})\,.$$
Other than $k=1$, we exclude odd $k$ because $m↓k=0$; henceforth, $k$~is
even. Then $g↓k(x)=g↓{k+1}(x)$, so
$$E↓k(a,b)=E↓{k+1}(a,b)=\int↓a↑bf↑{(k+1)}(x)q↓{k+1}(x)\,dx\,.$$

Suppose $f↑{(k+2)}$ and $f↑{(k+4)}$ have the same (constant) signs in $(a,b)$.
We shall see that $E↓k(a,b)$ has absolute value no greater than
$T↓{k+2}(a,b)$, and has the same sign; that is, the error is dominated
by the first omitted term.

Let $m$ be an even integer, $c$ any integer in $a≤c≤b-1$. Then
$$\eqalign{E↓m(c,c+1)&=E↓{m+1}(c,c+1)\cr
\noalign{\smallskip}
&=\int↓c↑{c+1}f↑{(m+1)}(x)q↓{m+1}(x)\,dx\cr
\noalign{\smallskip}
&=\int↓0↑1f↑{(m+1)}(c+x)q↓{m+1}(x)\,dx\cr
\noalign{\smallskip}
&=\int↓0↑1f↑{(m+1)}(c+1-x)q↓{m+1}(1-x)\,dx\cr
\noalign{\smallskip}
&=-\int↓0↑1f↑{(m+1)}(c+1-x)q↓{m+1}(x)\,dx\cr
\noalign{\smallskip}
&={1\over 2}\,\int↓0↑1\bigl(f↑{(m+1)}(c+x)-f↑{(m+1)}(c+1-x)\bigr)q↓{m+1}(x)\,dx\cr
\noalign{\smallskip}
&=\int↓{1/2}↑1\bigl(f↑{(m+1)}(c+x)-f↑{(m+1)}(c+1-x)\bigr)q↓{m+1}(x)\,dx\cr}$$
where the first factor has the sign of $f↑{(m+2)}(\zeta)$.

Substituting $m=k$ and $k+2$, and summing, $E↓k(a,b)$ has the same sign
as $-E↓{k+2}(a,b)$. Add the three equations
$$\displaylines{T↓{k+2}=(g↓{k+2})↓a↑b-(g↓k)↓a↑b\cr
\noalign{\smallskip}
\sigma↓{ab}=(g↓k)↓a↑b+E↓k(a,b)\cr
\noalign{\smallskip}
(g↓{k+2})↓a↑b+E↓{k+2}(a,b)=\sigma↓{ab}\cr}$$
to find $T↓{k+2}=E↓k(a,b)-E↓{k+2}(a,b)$. Therefore $E↓k(a,b)$ is a fraction
of~$T↓{k+2}$. This gives a powerful error bound:
{\it If ${\rm sign}(f↑{(k+2)})={\rm sign}(f↑{(k+4)})$ is constant in $(a,b)$,
the error in Euler's approximation is a fraction of the first omitted
term.}
For even $k$, it is known that $|m↓k|$ is a little under $3\times 6↑{-k}$.

Warning: Euler's series $g↓k(x)$ is inherently finite. Taken as an infinite
series, it almost always diverges to infinity. Typically, it is used
with $k=2,4$, or~6.

The bound above is useful for summing series like $x↑{-d}$ and $x↑{-d}\ln x$.

To calculate $\sigma↓{ab}=\sum↓{a<i≤b}f(i)$ for large $b$, with error
no more than~$\epsilon$, use this technique. Choose $k$ in Euler's
formula such that $\int↓a↑∞|f↑{(k-1)}(x)|\,dx$ is finite. As a practical
matter, choose $f↑{(k-1)}(x)\in O(x↑{-4})$. Next choose~$n$ just large
enough that $|E↓k(\alpha,\beta)|<\epsilon$ for all $\alpha,\beta ≥n$.
Then (Variant~1) $\sigma↓{ab}=\sigma↓{an}+\sigma↓{nb}=\sigma↓{an}-g(n)
+g(b)+E↓k(n,b)$. Calculate $\gamma↓n=\sigma↓{an}-g(n)$, noticing its
independence of~$b$; the result is $\sigma↓{ab}\approx \gamma↓n+g(b)$, with
error bound $|E↓k(n,b)|<\epsilon$.  (Variant~2) Obtain $\gamma↓∞=\lim↓{c→∞}
\bigl(\sigma↓{ac}-g(c)\bigr)$. 
Because of the choice of~$k$, this limit exists.
$\sigma↓{ab}=\sigma↓{ac}-g(c)+g(b)-E↓k(b,c)$, so, taking limits,
$\sigma↓{ab}\approx\gamma↓∞+g(b)$, with error bound $|\lim↓{c→∞}E↓k(b,c)|<\epsilon$.

\vfill\eject

\line{\bf Examples.\hfil}

$$\vcenter{\halign{$\hfil#\hfil$\quad
&$\hfil#\hfil$\quad
&$\hfil#\hfil$\quad
&$\hfil#\hfil$\cr
f(x)&1/x&\ln x&x\ln x\hbox{ (Knuth 1.2.11.2 Ex.~7)}\cr
\noalign{\smallskip}
f↑{(-1)}(x)&\ln x&x(\ln x-1)&x↑2\ln x/2-x↑2/4\cr
\noalign{\smallskip}
f'(x)&-1/x↑2&1/x&\ln x+1\cr
\noalign{\smallskip}
f'''(x)&-6/x↑4&2/x↑3&-1/x↑2\cr
\noalign{\smallskip}
f↑{(5)}(x)&-120/x↑6&24/x↑3&-6/x↑4\cr
\noalign{\smallskip}
f↑{(7)}(x)&-5040/x↑8&720/x↑7&-120/x↑6\cr
\noalign{\smallskip}
\epsilon&10↑{-10}&10↑{-10}&2\times 10↑{-8}\cr
\noalign{\smallskip}
k&4&6&4\cr
\noalign{\smallskip}
n&19&10&10\cr
\noalign{\smallskip}
g↓k(x)&\ln x+{1\over 2x}%
&x(\ln x-1)+{\ln x\over 2}+{1\over 12x}&\ln x\bigl({x↑2\over 2}+{x\over 2}%
+{1\over 12}\bigr)-{x↑2\over 4}\cr
\noalign{\smallskip}
&\null-{1\over 12x↑2}+{1\over 120x↑4}&\null-{1\over 360x↑3}+{1\over 1260x↑5}%
&\null+{1\over 12}+{1\over 720x↑2}-{1\over 5040x↑4}\cr
\noalign{\smallskip}
g↓k(n)&2.970523996&14.18547404&101.9174094\cr
\noalign{\smallskip}
a&0&1&1\cr
\noalign{\smallskip}
\sigma↓{an}&H↓{19}=3.547739657&\ln(10!)=15.10441257&2\ln 2+\cdots +10\ln 10\cr
\noalign{\smallskip}
&&&=102.0828306\cr
\noalign{\smallskip}
\gamma↓n&0.5772156651&0.91893853&0.16542119\cr
\noalign{\smallskip}
\gamma↓∞&0.5772156649&0.9189385332&\hbox{Glaisher's constant $-1/12$}\cr
\noalign{\smallskip}
&=\hbox{ Euler's constant}&=\ln\sqrt{2\pi}\cr
\noalign{\smallskip}
b&1000&50&100\cr
\noalign{\smallskip}
g↓k(b)&6.908255196&147.5588284&20756.57654\cr
\noalign{\smallskip}
\gamma↓n+g↓k(b)&7.485470861&148.4777669&20756.74196\cr
\noalign{\smallskip}
\gamma↓∞+g↓k(b)&7.485470861&148.4777669\cr
\noalign{\smallskip}
&=H↓{1000}&=\ln(50!)\cr
\noalign{\smallskip}
&&\hbox{so $50!=3.041409\times 10↑{64}$}&\hbox{so $2↑2\times 3↑3\times\cdots
\times 100↑{100}$}\cr
\noalign{\smallskip}
&&&=3.4554\times 10↑{9014}\cr
}}$$
In retrospect, since the rounded results were only good to 10SD, the $\epsilon$'s
should have been $10↑{-9}$, $10↑{-7}$, and $10↑{-5}$; the $n$'s could have
been 13, 4, and~3.

\vfill\eject

\line{\hfil Odd bits \# 1}

An alternating series can be treated as the difference of two series of
constant sign, to which Euler's formula is applied:
$$\eqalign{1-{1\over 2}+{1\over 3}&-{1\over 4}+\cdots -{1\over 2m}\cr
\noalign{\smallskip}
&=\left(1+{1\over 2}+{1\over 3}+\cdots +{1\over 2m}\right)-2\left(
{1\over 2}+{1\over 4}+\cdots +{1\over 2m}\right)\cr
\noalign{\smallskip}
&=H↓{2m}-H↓m\cr
\noalign{\smallskip}
&=\left(\ln(2m)+{1\over 4m}-{1\over 48m↑2}+{1\over 1920m↑4}\right)\cr
\noalign{\smallskip}
&\quad\null -\left(\ln m+{1\over 2m}-{1\over 12m↑2}+{1\over 120m↑4}\right)
+E↓4(m,2m)\cr
\noalign{\smallskip}
&=\ln 2-{1\over 4m}+{1\over 16m↑2}-{\theta\over 128m↑4}\,,\quad
0<\theta<1\,.\cr}$$

\bigskip
\bigskip
\line{\hfil Odd bits \# 2}

The formulas derived above may be applied mechanically to certain non-integer
instances. For example, treat ${1\over 1.1}+{1\over 2.1}+\cdots +{1\over 100.1}$
as $\sigma↓{0.1,100.1}$ for $f(x)=1/x$. Use $k=4$, $n=13.1$, $g↓k(x)$ as
before, $\sigma↓{0.1,13.1}={1\over 1.1}+{1\over 2.1}+\cdots +{1\over 13.1}$,
and calculate the result
$$\sigma↓{0.1,13.1}-g↓k(13.1)+g↓k(100.1)\,.$$

\vfill\eject

\line{[Needs work] [Redundant?]\hfil Odd bits \# 3}

Let
$$\eqalign{Q&=\left({2\cdot 4\cdot 6\cdot\ldots\cdot (2m)\over 1\cdot 3\cdot 5\cdot
\ldots\cdot (2m-1)}\right)\,,\cr
S&=\ln Q=-\ln 1+\ln 2-\ln 3+\cdots -\ln(2m-1)+\ln(2m)\cr
\noalign{\smallskip}
&=-(\ln 1+\ln 2+\cdots +\ln 2m)+2\bigl(\ln 2+\ln 4+\cdots +\ln(2m)\bigr)\cr
\noalign{\smallskip}
&=-(\ln 1+\cdots +\ln 2m)+2m\ln 2+2(\ln 1+\cdots +\ln m)\cr
\noalign{\smallskip}
&=-\bigl(\gamma↓∞+g↓2(2m)\bigr)+2m\ln 2+2\bigl(\gamma↓∞+g↓2(m)\bigr)
+O(1/m↑3)\cr
\noalign{\smallskip}
&=\gamma↓∞-2m(\ln 2m-1)-{\ln(2m)\over 2}-{1\over 24m}\cr
\noalign{\smallskip}
&\quad\null+2m\ln 2+2m(\ln m-1)+{2\ln m\over 2}+{1\over 6m}+O(?)\cr
\noalign{\smallskip}
&={\ln 2\pi -\ln 2\over 2}+{\ln m\over 2}+{1\over 8m}+O(?)\cr
\noalign{\smallskip}
&={\ln (m\pi)\over 2}+{1\over 8m}+O(?)\cr
\noalign{\smallskip}
Q&=\exp(S)=\sqrt{m\pi}\cdot\left(1+{1\over 8m}+{1\over 128m↑2}\cdots\right)
+O(?)\cr
\noalign{\smallskip}
&=\sqrt{m\pi}+O(m↑{-1/2})\,.\cr}$$

{\rmn
\disleft 38pt::
{\bf Exercise.}
Show $Q↑2/m=\pi+O(1/m)$. Show by regrouping that
$${8\over 9}\cdot {24\over 25}\cdot {48\over 49}\cdots ={\pi\over 4}\,.$$
(The above assumes $\gamma↓∞=\ln\sqrt{2\pi}$.)

}





\bigskip
\parindent0pt
\copyright 1987 Robert W. Floyd.
First draft (not published)
May 21, 1987.
%revised date;
%subsequently revised.

\bye